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LARGE EOOT SIMULATION OF TURBULENT CHANNEL FLOW - 
ILLIAC IV CALCULATION 

John Mm* and Parvfz Moln* 

Ames Research Center, NASA, Moffett Field, California 94035, U.S.A. 


SUMMARV 

The three-dimensional time-dependent equations of motion have been numerically Integrated for fully- 
developed turbulent channel flow. The large-scale flow field Is obtained directly from the solution of 
these equations, and the small-scale field motions are simulated through an eddy viscosity model. The ealeu 
lotions are carried out on the ILLIAC IV computer with 64 * 64 * 64 grid points. 

The computed flow patterns show that the wall layer consists of coherent structures of low-speed and 
high-speed streaks alternating in the spanwlso direction. These structures were absent In the regions away 
from the wall. Hot spots, small localized regions of very large turbulent shear stress, are frequently 
observed. Very close to the wall , these hot spots are associated with u" >0 and v < 0 (sweep) i away from 
the wall, they are due to u" < 0 and v > 0 (burst). The profiles of the pressure velocity-gradient correla 
tlons show a significant transfer of energy from the normal to the spanwlso component of turbulent kinetic 
energy in the Immediate neighborhood of the wall ("the splattlng effect"). 


NOMENCLATURE 



The overbar (“) denotes the filtered component and the 

prime (•) denotes subgrid scale (SGS) component, 

c s 

Smagorlnsky's constant 

u" 

2 u - <u> 

G(x 

- x 1 ) filter function 

u i 

velocity In the 1 -di roctlon 

h 1 

mesh size in the 1 -direction 

if 

Fourier transform of u^ 

V 

. Vr 

V 

U T 

shear velocity « t'TJo 

velocity in the vortical direction 

k 

wave numbc" e /k]-’ + k 3 z 

V 

k 1 

wave number In the 1-dtrectlon 

W 

velocity in the spanwlse direction 


X, Xj 

streamwlse coordinate 

L 

length nf the computational box In the 

coordinate In the 1-dlrectlon 

X 

x-dlrectlon 

*1 

L z 

length of the computational box in the 
z-dlrectlon 

X, *' 

coordinate vector 

i 

SGS length scale 

y- x 2 

coordinate In the direction normal to the walls 

N 

number of mesh points In the y-dlrectlon 

*W 

distance to the nearest wall 

P 

pressure 

y + 

Wt 

V 

P 

. P . R kk 
' o + T 

z. x 3 

spanwlse coordinate 


P , 1 z-rr- , R kk 

c 1Jk 

the completely antisymmetric tensor of rank 3 

P* 

3 P + I u j u j + T 

X 

mean streak spacing 

P 

Fourier transform of p 

*1 

mean spacing of the turbulent structures In the 
1-dlrection 

<1 

root-mean-square velocity 

x, * 

. 

Re 

Reynolds number based on channel half- 

*1 

5 V 


width and the centerline velocity 


Xu 


*+ 

- T 

Re 

Reynolds number based on channel half- 


V 

T 

width and shear velocity 


jth meshpoint in the vortical direction of the 


.... . 

transformed (uniform mesh) space 

R 1j 

3 VV + u j’ 5 i *Yi' 

p 

density 

S 1J 

, /au . Du A 

a 1 (axj + axT| strafn r3te tensor 

T 1j 

n WlJ 

R 1j ' 3 



T w 

mean wall sherr stress 

t 

dimensionless time 




A t 

dimensionless time step 

u 

streamwlse velocity 
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V 

kinematic viscosity 

< > horizontal average 

V T 

eddy viscosity 

< >*■ time average 

"i 

vortfetty in the 1-direction 

Subscripts 

B x 

vortfuity in the x-direction 

w wall value 


J 1 H 

SGS subgrid scale 

a ij 

(o w 

Superscript 



n time step 

1. 

INTRODUCTION 



The technique of 1 ai*ge eddy s tniuTatf on (LES) Is a relatively new method for computing turbulent flows, 
ine primary motivation for its undertaking is that the large eddy turbulence structures are clearly flow- 
dependent (e.g,, Jets vs boundary layers) and hence they are difficult if not impossible to model. On the 
other hand, there is experimental evidence {e.g., Ref. )) that small eddies are universal in character, and 
consequently much more amenable to general modeling. 


In LES, the large-scale motions are computed directly using three-dimensional time-dependent computa- 
tion, and the small-*'’*'” ■">*■<»“« n._ a .. — , 

derived by averaging 


in, and the small-scale motions are modeled. The dynamical equations for the large-scalf field are 

, - ^raging the Navier-Stokes equations over volumes in space that are small compared to the ovoral 
dimensions of the flow field. This averaging is to provide sufficient smoothing of the flow variables, so 


they can be represented on a relatively coarse mesh. The resulting equations for the largo eddies contain 
terms that involve small-scale turbulence. These terms are replaced by models that are to represent the 
interaction between the resolved and unresolved (subgrid scale, SGS) field motions. 


One of the most extensive applications of LES has been to the problem of decay of homogeneous isotropic 
turbulence (see Refs. 2-4). A variety of numerical methods and subgrid-scale turbulence models was incorpo- 
rated to compute this flow. Both the pressure-velocity and the vorticity-stream function formulations of 
the dynamical equations were used. These studies have shown that homogeneous turbulent flows can be reason- 
ably simulated using simple eddy-viscosity models. 


The first application of LES was made by Oeardorff (Ref. 5), who simulated a fully developed turbulent 
channel flow at a very large Reynolds number. Utilizing a modest number of grid points (6,720), he showed 
that three-dimensional numerical simulation of turbulence (at least for simple flows) Is feasible. His 
calculations predicted some of the features of turbulent channel flow with reasonable success and demon- 
strated the potential of LES for prediction and analysis of turbulent flows. 


Schumann (Ref, 6) has also performed numerical simulation .f turbulent channel flow. In addition, lie 
has applied LES to cylindrical geometries (annuli). Ha used up to 10 times more grid points than Oeardorff 
and a much more complex subgrid-scale model. In that model, an additional equation for SGS turbulent kinetic 
energy was Integrated, However, the results showed no significant improvement over the case in which eddy- 
viscosity models wore used (Ref. 6). 


In the calculations of channel flow described above, no attempt was made to compute the flow In the 
vicinity of the walls, A great portion of turbulent kinetic energy production takes place in this region 
(soo Ref. 7). Therefore, by using artificial velocity boundary conditions well beyond the viscous sublayer 
and buffer layer, a significant fraction of the dynamics of turbulence in the entire flow was effectively 
modeled. In addition, it should be noted that the boundary conditions used in the latter calculations 
assume that in thq log layer, the velocity fluctuations are in phase with the wall shear stress fluctuations. 
This assumption is not supported by experimental measurements (Ref. 8). 

Moin et al. (Ref, g) simulated the channel flow, Including the viscous region near the wall. The exact 
no-slip boundary conditions were used at the walls. In their computations, only 16 uniformly spaced grid 
points were used in each of the streamwise (x) and spanwise (zj directions and 65 nonuniformly spaced mesh 
points were used In the y-directlnn. The grid resolution was especially Inadequate In the z-direction to 
resolve the now well-known streaky structures In the vlcinfty of the wall. In spite of this, computations 
old display some of the well-established features of the wall region. In particular, the results showed 
coherent structures of low-speed and high-speed fluid alternating ii the viscous region near the wall, though 
not at their proper scale. The overall agreement of the computed mean-velocity profile and turbulent statis- 
tics with experimental data was satisfactory. 


Encouraged by the results of the above coarse calculation, the present numerical simulation of channel 
flow with 262,144 grid points (64 x 64 « 64) was undertaken. The ILLIAC IV computer, a parallel processor, 
was chosen for this purpose. Although the grid resolution In the spanwise direction is still not sufficient 
for an adequate representation of the wall-layer streaks, it is a significant improvement over the earlier 
calculation. This, in turn, allows a more realistic and comprehensive study of the structure and mechanics 
of this flow. 


This paper Is the result of a work that is now In progress and is essentially Intended to demonstrate 
some of the capabilities of LES in the prediction and analyses of wall-bounded turbulent shear flows. In 
Sec. 2, the dynamical equations for large-scale field motions are derived. The subgrid model that was used 
is described in Sec. 3; Section 4 describes the computational grid network and its relation to the observed 
physical length scales in the flow. The numerical method Is briefly outlined in Sec, 5; the data management 
process is taken up in Sec. 6; and in Sec. 7, we examine some aspects of the mechanics and structure of the 
flow, both in the vicinity of the wall and in regions away from the wall, and an attempt is made to correlate 
numerical results with laboratory observations. In Sec. 8, we present the computed flow statistics, which 
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include the mean-velocity profile, turbulent intensities, and turbulence shear stress. In that section, we 
will point out some of the deficiencies of the subgrid-scale model used and suggest improvements. Finally, 
conclusions are presented In Sec. 9, 

Z. GOVERNING EQUATIONS FOR THE LARGE-SCALE FIELD 

Tne first step in LES is the definition of the large-scale field. Each flow variable f is decomposed 
as follows! 


f - f + f' 


(1) 


Here, the overbar denotes the large-scale or "filtered" field and the prime indicates the residual or "sub- 
grid" field. Following Leonard (Ref, 10) wo define the large-scale field as: 

( 2 ) 


f(x) « J* G(x,i')f(x.')dx' 


where G is the filter function and the integral is extended over the whole flow field. In the horizontal 
ptanes (x-z), several possible choices for the filter function are available. Unless otherwise stated, most 
of the calculations reported here were carried out using a Gaussian filter, G(x-x‘ ,z-z'). The width of the 
Gaussian function characterizes the smallest scales of motion retained in the filtered field (the largest 
scales in the residual field). He assume that the filtering in the planes parallel to the walls provides 
sufficient smoothing in the vertical directions as wollj our computations support this assumption. In addi- 
tion, it should be noted that we use second-order finite difference schemes to approximate partial derivatives 
In the Xj-direction and such schemes have an implicit filtering effect associated with them. For further 
details see Hofn ct al. (Ref, 9). 

After applying the filtering operation (Eq, (Z)) to the incompressible Navfer-Stokes and the continuity 
equations, the governing equations for the filtered field may be written 


^Uj 

3T 


c ijk y 


auj 

3x7 


. isl + 


»X 


» 0 


3*u, 


ii ' 3Xj T iJ + ReT axj3Xj 


(3) 


(4) 


where we have decompo^d as in (1) and 




pqk ax 


T U = R 1J 


. R kkf» 


R,, p u.'u,' + u.'ui + u.u 


'1J U 1 U J 


j 


yi 


5 1 r-r— + R|fl* _ _ _ 

P + ? u j u J TT 3 p + 1 u j u j 




Here, the variables are nondfmensional using the channel half-width 6 and the shear velocity u T ,'r^/p. 
The calculations will be carried out for a fixed streamwise mean-pressure gradient which is accounted for 
by the 4fj term in the momentum Eq. (3). 

3. RESIDUAL STRESS MODEL 

The remaining unknown quantity in Eq. (3) is t-h. This tone represents the subgrid-scale stresses and 
must be modeled. In the present calculations we have used an eddy viscosity model, 


T ij “ - Zv T S 1j 


where 


a a 


l 

1 


3Xj 


(5a) 


(Sb) 


The small-scale eddy viscosity vt represents the action of the unresolved scales of motion on the 
resolved scales. Hence, as the resolution gets better, should get smaller. This suggests that vj 
should scale on a length scale t which is directly related to the computational resolution. The model 
most commonly used for v T and the one we use here is the Smagorinsky model, 


v T - M’-'TC 


( 6 ) 


where C s » 0.1 (Ref, 5) is a dimensionless constant and i is a dimensionless representative of the grid 
resolution, here assumed to be (Ref. 5): 

t r (hj • ha(y) • tu) 1 / 3 


(7) 


14-4 


This expression for t Is probably appropriate 
anisotropy (Ref, 6), In the present calculation 
the vicinity of the walls, and hence use of Eq, 

Insight into the role of i and to help guide 1 
with a modification described below. 

Near the walls, the subgrid-scale turbulence Reynolds number, defined as 



1 s very sma 1 1 , and hence one expects the value of the eddy viscosity coefficient to bo very small, In our 
calculations, we have found that the damping provided by the presence of (hjly)) 1 ' 3 In Eq. (7) Is not suffi- 
clent, and excessively large subgrid-scale stresses are formed near the wall. Therefore, In. the present 
calculations we have multiplied l (Eq. (7)) by an exponential damping function 1 - exp(-y + /50). 

The eddy-viscosity model used here is best rationalized for isotropic turbulence at the scale of the 
computational grid. The fundamental assumption behind this model Is that the resolution scale lies within 
an inertia) range with the -5/3 power spectrum (Ref. 11). It is clear that for the moderate Reynolds number 
,1 j ' if 1 ? c V e are considering and the nature of the grid volumes used, the above assumptions are not 
satisfied. This is particularly true in the highly viscous region in the vicinity of thn walls. Thus, the 
present simulation is Viewed as a challenge to the eddy-viscosi ty model used. 

A critical test for the largo eddy simulation technique is the prediction of the logarithmic layer 
and the von Karman constant. 1 This is one of the reasons for not utTI7zTng~ti}io mixing- length model in the 
present calculations to account for inhomogenef ty due to the moan shear (Ref. 6). Such a model is known to 
postdlct" the correct mean-velocity profile. 

4. THE COMPUTATIONAL GRID 

The availability of computer resources restricts the size of calculations possible. For a given number 
of grid points N, wo have to choose the grid size(s) based on the known physical properties of turbulent 
channel flow under consideration. 

In the vertical direction (-1 < y 1), a nonuniforra grid spacing is used. The following transformation 
gives the location of grid points in the vertical direction (Ref. 9): 

yj * j tanh [Cj tanh'UaJ] ( 9 ) 

where 

Cj ■ 1 + Z(J - 1 )/(N - 2) (10) 


J « 1,2 N 

N Is the total number of grid points in the y direction, and the adjustable parameter of transformation is a 
(0 < a < 1). We used a “ 0.98346, N » 64. This value of a was selected so that the above grid distribu- 
tion in the y-direction is sufficient to resolve the viscous sublayer (y + < 5). 

The length L x and Lz of the computational box in the stroamwlso (x) and spanwisc (z) direction, in 
which periodic boundary conditions are used, should be long enough to include the important large eddies 
(Refs. 6, 12). Based on the two-point correlation measurements of Comte-Bellot (Ref. 13). we used Lx ■ 2n, 
a , *-z “ 4»/3. We have used 64 uniformly spaced grid points in each of the streamwiso and spanwisc direc- 
tions. With the above choices for L x and Lz. the nondimonsional grid spacings in the horizontal directions 
expressed in the wail units are; 

h l + » 63 
hj + « 42 

In the wall region, the importa it large eddfes are the "streaks” (Ref, 14). These have a moan spanwisc 
spacing corresponding to => 100. It is clear that our grid resolution in the spanwise direction is not 
quite sufficient to resolve the stro/.ks. This is especially true when we note that the above value for *, + 

Is based on on ensemble of measurements, and at a given instant streaks with a finer spacing than * 3 + can 
be formed. As we shall see, however, calculations did reveal these structures, though not at their proper 
scale. 

With relatively minor modifications to the present computer program, we are able to perform calculations 
with 64 x 64 x 128 grid points in the x, y, and z directions, respectively. It is expected that in this 
simulation the spacing of the wall-layer streaks will be more In line with the laboratory observations. 

5. NUMERICAL METHOD 

A complete description of the numerical method used Is given in Ref. 15. Here, we give a brief outline 
of the method and minor modifications that were made to enhance the data management process. The partial 
derivatives In the x 2 direction were approximated by second-order central difference formulae. In the 
X] and Xj directions, partial derivatives were evaluated pseudospoctrally (Ref. 16). With a given number of 
grid points, the use of the pseudospoctral method in any given direction gives us the best possible resolution 
in that direction. This Is particularly useful in the x 3 direction whore wo face a lack of grid resolution 
(Sec. 4), 


only for cases in which there is no significant grid 
, tho computational grid is very elongated (hi.h, » h 2 ) in 
(7) is not strictly justified. However, to gain a better 
ts selection in future calculations, wo have used Eq. (7) 
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Tine advancement is made using a seml-lmpltcit method, Pressure, viscous tonus, and part of the subgrid- 
seale model are treated Implicitly, whereas explicit time advancement Is used for the remaining nonlinear 
terms. The equation of continuity Is solved directly, Second-order Adams Basnfortli (Ref. 17) and 
Crank-Nicolson (Ref. 10) methods are used for explicit and Implicit time advancement, respectively. 

Next, we Fourier transform the resulting equations In Xj and Xj directions. This converts the above 
sot of partial differential equations to the following sot of ordinary differential equations for tliu variables 
at tin® stop n + 1, for every pair of Fourier wave numbers k 3 and k 3 , with y » x 2 as the Independent 
variable, 


»*U? H + 
■ jy* 

(D, - 

k J )ui +l ♦ 

1k|K| f 

p" + ‘ 

“ Qi” 

Ola) 

a-’uf 
’ ay” 

+ (P; 

- k^Of 1 

.. At 8^ 

♦ Pa r ay” 

- Qz" 

Olb) 








3 Uj 

(#3 * 

war + 

•Ml / £- 

p n+1 

- Qi n 

(lie) 

ay* 










au," M 
’ ay + 


ik 3 ar 


(lid) 


Here, t.| (1 * 1,2,3) aro known functions of Ro T and s^, and s n represent the terms Involving the 
velocity and pressure field at time-step n and n - 1 (see Ref, 15; 


In addition to the use of Implicit time advancement for all the viscous terms, the algorithm used In the 
present study Is different In one other respect from the one ‘described In Ref. 15. For reasons that will be 
explained In tne next sec tier., Eqs. (11a) and (11c) were mul tf pi led by lk t and Ik*, respectively. Thus, 
the dependent variables for the time-advancement equations are Ikju, v, and lk_,w rather than u, v, and w. 

The remaining stops In the solution procedure are ns follows. Finite difference operators (described 
above) are used to approximate 3/ay and a*/3y 2 . This gives a sot of linear algebraic equations for the 
Fourier transform of dependent variables. This system Is of block tridiagonal form and can be solved very 
efficiently. No-slip boundary conditions are used at the solid boundaries. Finally, inversion of the 
Fourier transform gives the velocity and pressure field at time-step n + l. 


The initial velocity field was the same as the one used In Ref, 9 Interpolated on the finer grid used 

here. 


6, DATA MANAGEMENT 

In large simulations, the high-speed random-access memory of the computer on hand may not hold the 
entire data base of the problem being considered. In the present case, the core memory of the 1LUAC IV Is 
large enough to hold only a few planes of velocity pressure field. Therefore, It Is essential to manage the 
flow of data efficiently between the core memory and the disk memory where the entire data base resides. In 
general, separate passes over the data base are requl 'ed for each time step and the task Is to minimize the 
required number of such passes. The following describes a data management process employed In the present 
simulation. 

The system of Eq, (11) must he solved for both real and Imaginary parts of the dependent variables. 

This necessarily means that two_passes through the data base are required: one for real parts of u 3 and fit 

and Imaginary parts of u 2 and p, and the other for Imaginary parts of iij and u 3 and real parts of u 2 and p. 

To avoid an extra pass through the data base, we multiply Eqs, (11a) and (11c) by ikj and lk 3 , respec- 

tively (Ref. 19). (These multiplications In Fourier space amount to differentiations In real space.) 


,2,-jin i 

‘ 1 + (B, - k*)uT - - Q,» 

ay* 

02a) 

^MPz-k ^ 1 .q 2 n 

ay- J 

02b) 

sZii n ' H 

' 3 + (P 3 - k 2 )u 3 +l - k 3 2 B 3 4rp ,,+1 » Q, 11 

ay 2 

02c) 

,.-nH 

02d) 

ik 3 (l 3 i Qi" » lk ; Q, n i i) 2 n = Q 2 "; and q 3 " ■ lk 3 q 3 ". 

The above system of 
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equations can be solved with one pass through the data base, but two extra integrations fn the Fourier space 
are required to obtain u, and u 3 in physical space. It should be noted, however, that such integrations cost 
far less than an I/O /ass. In addition, to avoid the loss of Information, upon differentiation, the Fourier 
mode associated with a null wave number is simply not multiplied hy its wave number (i.e., zero) and, 
similarly, it is not divided by its wave number upon integration. This implies that ut, 0?, and u« in 
Eqs, 12 should be understood as 


ujlO.y.ks); ikiUi(k J ,y,ifj),k I / 0 

uj(ki,y,k 3 ) 


iij(fcj.y.O)t ikjiijfkj.y.kjJ.kj / 0 

The system of Eqs. 02) is. solved by two separate passes through the data base, In PASS 1, the right* 
hand sides of these equations, Qj (1 ■ 1,2,3), are evaluated and in PASS 2, the block trldlagonal system is 
solved. To compute the right-hand side vector in PASS 1, differentiations in all spatial directions are 
required. Since the pseudospectral method Is used in the horizontal directions (x and z) and a finite- 

difference scheme is used in the normal direction {central difference), all the data in an (x - z) plane are 

needed for operators in these directions and the data for at least three adjacent pianos are needed for finite 

difference operators, in the y direction. Therefore, in PASS 1, two (x - z) planes are brought into the core 

to be handled by a double buffer scheme. One complete pass through the data base is required to complete 
PASS 1 , 

In PASS 2, the block tridiagonal system must be solved for each k 3 and k 3 . In this pass, two (y - k 3 ) 
planes are brought into the core, A special algorithm had to be> developed to solve the block tridiagonal 
matrix because of the limitation on the core size. In a conventional -block tridiagonal solver, all the 
results of forward sweep are stored to be used in backward sweep. For the present simulation, this would 
require half of the total core size {i.e., 16 * 64 x 64) which is not feasible. Hence, a special algorithm* 
was developed so that only a part of the results of the forward sweep is stored in the memory and the rest is 
recomputed as necessary in the backward sweep. Although this requires extra computations In the backward 
sweep, this method is much more efficient than performing the extra 1/0 passes that would otherwise have 
been necessary. 

The computation described here was carried out on the ILLIAC IV computer at Ames Research Center. The 
dimensionless time step, during most of the calculations, was set at at « 0.001, The computer time per 
time-step (CPU and I/O time) was about 22 sec. This computational speed has been achieved with a full use 
of the parallel processing capabilities of the ILLIAC IV and the data management process Just described, 

7. DETAILED FLOW STRUCTURES 

In this section, we investigate the detailed flow patterns by examining contour plots of typical 
Instantaneous velocity and vorticity fields in x-z, x-y, and y-z planes. In all these plots positive 
values are contoured by solid lines and negative values are contoured by dashed lines. In addition, all 
tlte plots are obtained at a given dimensionless time (t = 1,4). 

Figure 1 shows patterns of u" in an x-z plane very close to the lower wall (y + - 16.1). The striking 

feature of this figure is the existence of highly elongated (in the x-ulrection) •eglons of high-speed 

fluid located adjacent to low-speed ones. This picture of the flow pattern in the vicinity of the wall is 
in agreement with experimental observations (Refs. 20, 21) that the wr.ll layer consists of relatively coherent 
structures of low-speed and high-speed streaks alternating in the spanwise direction. Examination of the 
typical spanwise spacing of these structures shows significant Improvement over the earlier simulation 
(Ref. 9) where only 16 uniform grid points were used In each of the spanwise and streamwlse directions. How- 
ever, the typical spacing of these streaks is still about 3 times larger than the experimentally observed 
moan value of x 3 + » 100. This is expected, since our computational grid size In the spanwise direction is 
too largo to resolve the wall layer streaks in their proper scale {Sec, 4). 

Figure,^ shows patterns of u" in an x-z plane far away from the wall (y/i » 0.73), It is clear 

that the u patterns in the regions away from the wall do not show the coherent streaky structures that are 

characteristic of wall-layer turbulence. This is also In agreement with the experimental observations 
(Rof. 20), In fact, it is difficult to associate a definite structural pattern to u in the regions away 
from the wall. 

Since turbulent energy production is directly proportional to -<uv> t , it is Important to study the 
instantaneous map of u"v, Figure 3 shows the patterns of u"v In the same x-z plane as in Fig. 1; 
that is, very close to the wall (y 1 *' * 16,1). Examination of this figure reveals several points related to 
the dynamics of wall-layer turbulence that deserve attention. First, it can bo seen that the regions with 
negative u"v, which have a positive contribution to the production of average turbulent kinetic energy, 
constitute the overwhelming majority of the entire plane. Second, pronounced stroamwise elongation, the 
characteristic of the wall layer u“ eddies, is absent in 0"v patterns. This indicates that In contrast 
to u" eddies, v eddies are not significantly elongated in the x-direction. Third, there are several 
small regions (hot spots), that are associated with very large values (large concentrations of dashed lines) 
of -D"v. These regions are highly localized In spac-:, Overlaying Fig. 3 on Fig, 1 reveals that the great 
majority of the "hot spots" are associated with u" > 0 (hence, 5 « 0). Thus, It appears that In the close 
vicinity of the wall most of the regions with very large values of (-u"v) are associated with hif.h-speed 
fluid approaching the wall (sweeps) rather than low-speed fluid being ejected from the wall (bursts). With 
combined visual and hot-wire measurements, Falco (Ref. 22) has identified a new flow module in the vicinity 
of the wall. These relatively small but energetic structures (called pockets) appear to be footprints of 
high-speed fluid moving toward the wall. It Is possible that the hot spots identified here may be related 


*The original concept was suggested to us by Marshall Merriam, Ames Research Center. 
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to pockets. Figure 4 shows the contour plots of u v in tho x-z plane located at y + ■> 90, Examination 
of this figure and the corresponding u plot (not shown here) shows that in contrast tp the ncar-watl 
region most of the hot spots that can be identified In this plane are associated witii O' < 0 and v > 0, 
that is, with bursts. With quadrant analysis of uv, Brodkoy ct al. {Hof. 23) have found that most of tho 
contribution to -<uv>t in the wall region comes from sweeps, and that in the regions away from the wall it 
comes from ejections. This Is consistent with what is observed here in relation to Figs. 3 and 4. There 
are two other features in Fig, 4 that deserve attention. First, similar to Fig. 1, tho regions with negative 
u v constitute the overwhelming majority of the entire plane, Although there are regions with very largo 
positive u v, they are highly localized in spaco, Second, the maximum vah:» of (-B"v) in this.plane is 
17.81. This is about 20 times the expected <*uv»t at this plane. Such large excursions of u”v from Its 
expected mean value have been a frequent observation in tho laboratory (o.g., see Ref. 24). 

Figure 5 shows contour plots of u"v in an x-z piano far away from tho lower wall (y/a 0 0 ; 73). In 
contrast to planes located close to the lower wall (Figs, 3, a), where tho regions with negative u v domi- 
nated the entire planes, a significant portion of this plane Is associated with large positive u'v as well 
as negative u‘’v. The regions with the largest positive ii"v are associated with high-speed fluid moving 
toward the upper wall, ana tho regions with the largest -u"v seem to be evenly distributed among high-speed 
fluid movHgg toward the lower wall or low-speed fluid moving away from the lower wall. Finally, examination 
of the u v patterns in tho midplane (not shown hero) reveals that in contrast to the plane Just described 
(y/i « 0.73), the regions with the largest u“v are associated with bursts originating In the upper half of 
the channel, whereas the regions with the largast -u"v correspond to bursts originating in the lower half 
of tho channel, 

Among tho conceptual models of the inner region of turbulent boundary layers Is the streamwise vprtlclty 
model, fhls model portrays tho Inner region as being composed of pairs of long counter-rotating streamwise 
vortices located adjacent to each other, These long vortical structures, In turn, create low-speed and high- 
speed streaks alternating in the spanwise direction. Figure 6 shows the streamwise vorticity patterns in 
the same x-z plane as in Fig. 1 (y + ° 16). These patterns do not slnw elongated regions of positive and 
negative 3 X alternating in the spanwise direction. Moreover, no def' litc relationship appears to exist 
between the streak patterns shown in Fig. 1 and Z x patterns shown in Fig. 6. Therefore, the present simu- 
lation tends to dispute the validity of the vorticity model. 

Figures 7 and 0 show patterns of u" and in an x-y plane, z 0 15hj. For clarity, we have expanded 
the region 0 < y/a < 0.5, A pronounced feature of Fig. 7 is tho two regions of high-speed fluid (with res- 
pect to tho local mean velocity) that are inclined at oblique angles with respect to the wall. These struc- 
tures arc apparently associated with intense shear layers that are also inclined with respect to the wall 
(Fig, 81. Similar large-scale structures have also been observed in tho laboratory. Ft-om measurements of 
space-time correlation of wall shear stress and velocity fluctuations in a turbulent duct flow, Rajagopalon 
and Antonia (Ref. 8) have identified large-scale structures that are inclined at a moan angle of about 13° 
to the wall. At this time, we have not scanned a sufficient number of x-y pianos at widely spacod times 
to obtain the mean Inclination angle of these structures. 

in Figs. 9 through 14, contour plots of tho velocities and the streamwise vorticity in a y-z plane 
(x » 0) are shown. The contour plots in this plane reveal the existence of surprisingly well -organized 
structures in the wall region. Figure 9 shows a contour plot of f‘-,o streamwise velocity u". Note that the 
figure is stretched 4 times in the vertical direction and that the contour line patterns are thus distorted 
in that direction. Two important features can be observed in tills figure. First, away from the wall - for 
example, y/4 > 0.4 - no definite structure is discernible. Hear the wall, however, an alternating array of 
low-speed and high-speed fluid is noticeable. This array has a long streaky structure in the stroamwlse 
direction, as was shown in Fig. 1. Second, as we approach the wall, the size of the eddies decrease, 
gradually. Figure 10 is a magnified version of Fig. 9 close to tho wall, 0 ■; y + < 46. Again, the figure is 
highly stretched in the y direction so that the shapes of the flow structures are distorted. The array of 
low-speed and high-speed fluid is clearly discernible in this figure. This strikingly wel 1 -organ! zed flow 
structure In the wall region is consistent with the previous experimental observations (Ref. 20), although 
tho typical spacing between the streaks is no*, correct because of the insufficient spanwise grid spacings 
mentioned earlier. In addition to the well-oi eanized structure in the wall region, there exists a very 
intense shear layer in the vertical plane where the low-speed and high-speed fluids come close together. 

This could cause free-shear- layer- type Instabilities in this planes such instabilities might bo related to 
the experimental observations that the lifted stre.'ks oscillate not only in the vertical direction but also 
in the horizontal planes, 

Figure 11 shows a contour plot of tho normal veluc'ty v in the same plane as In Fig. 10. Here, a 
positive v (the solid lines) represents fluid moving «.>ny from the wall, and a negative v (the dashed 
lines) represents fluid moving toviard the wall. In tills figure we notice an array of fluid moving away and 
toward the wall. If we align Fig. 10 with Fig. 11, wo no-i-c that, generally, there exists a negative corre- 
lation between u“ and v. Note that in the vicinity of the ,'all, tho low-speed fluid elements (u < 0) are 
generally being ejected away from the wall (v > 0), while hi v v-spoed fluid elements are moving toward the 
wall, Clearly, tho fluid motions just described have a positive contribution to tho production of averaged 
turbulent kinetic energy. 

Figure 12 shows a contour plot of the opanwlv* velocity w, A positive w (solid lines) represents 
fluid moving to the right and a negative w (dashed line) represents fluid moving to the left. Note also 
that a significantly largo spanwise velocity gradient in y - that Is, ;.w/ay - exists due to the no-slip 
boundary conditions at tho wall. This results in substantial streamwise vorticity near tho wall, although 
flow Is not actually revolving in this region. We will come back to this later, if wo now align the con- 
tour plot of w with that of v, we can identify a definite flow pattern that exists In the wall region. 

A schematic illustration of this flow pattern is given in Fig. 15. This simplified illustration shows how 
low-speed streaks are being formed and lifted away from the wall. It is interesting to note that the rota- 
tion of the streamwise vorticity is in the opposite direction to the conventional vorticity model (Ref. 25) 
(see also Fig. 15b). 
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Figure 13 shows a contour plot of « x In the y-z plane at x * 0. It can be sen that S« Is 
concentrated only In the wall region, Away from the wall, tha strength of the vortldty becomes very weak 
and no organized structure Is discernible. Hear the wall, highly localized concentrations of S* appear, 
sometimes In pairs of opposite sign. Figure 14 Is a close-up of the wall region for > + < 46. Again, the 
figure Is highly stretched In^the vertical direction so that the pattern* are distorted. By comparing those 
contour pit vs with those of v and w, we can distinguish the streamwlse vortlclty associated with the revolv- 
ing fluid motion from the one associated with the velocity gradients. Recall that the cxiitencc of w* does 
not guarantee large-scale revolving fluid motion. In (act, most 5, very close tu the wait , say y f « 10, 
is due to 5w/oy and 1$ not related to the revolving action. Some of « x away from the wall, however, (e.g,, 
n the center In Fig, 14) Is associated with a large-scale reviving motlen. This Is In tgreeisent 
with the experimental observations. by flow visualization techniques (kef. 7) whom strong revolving motions 
ar<t observed away from the wall (y T > 10) and not very close to It. It should also be noted that although 
the strong vortical revolving fluid motion appears outside the sublayer, In the present simulation, the root* 
mean-square value of 5 X , <;,„*> 1/2 a | ways attains Its maximum at tha wall [ncto that ■ (Jw/ay)[waI1J. 

8. MEAN VELOCITY PROFILE AND TURBULENCE STATISTICS 

Figure 16 shows the mean-velocity profile <u> that has developed after two dimensionless time units. 

(One nondlmenslonal time unit corresponds approximately to the time In which a particle movlnq with center- 
line velocity travels 22a.) Mote that In the present study horizontal-average values are approximately 
ergodvc, The calculated velocity profile shows a distinct logarithmic region over an appreciable portion of 
the channel width. For comparison, we have also Included some of the available experimental data in this 
figure. The agreement of the computed mean-velocity profile with experimental data In mast of the channel 
Is satisfactory. In the vicinity of the wall, however, the values of the computed mean-velocity profile are 
rather low. This Is due to the presence of an excessively large eddy viscosity coefficient near the wall, 

To verify this observation, we carried out a set of calculations (starting from t * 1.0) where instead of 
the eddy viscosity model, we used a subgrid scale model similar to the one used by Fornborg (Ref, 26j 
in our numerical experiment, small-scale turbulence Is removed by a sharp cu* f' ; Uer at each time step). 
Although this model Is rather inadequate for proper representation of the Interact! 'i between the subgrid- 
scale and resolvable scale motions, ft suffices for our present purpose, especially ft the total time of 
Integration Is not large. Figure 17 (hows the resulting <u> profile at t * 1,5. It is clear that the 
profile of <u> has attained the proper values In the vicinity of the wall. In addition, the logarithmic 
*;y| r i U on £? analn e YldDnt. Figure 18 shows the profiles of resolvable normal turbulent Intensities. 

<u"*>i z, md <w J >i< i at the same time as In. Fig. 16, It can be seen that in agreement with 

oxpe rental measurements, generally, <u' l2 > t22 > <w 2 »i /* > <g*»i/2 throughout the channel, In addition, 

<u' >' 2 and <w‘> 1/2 attain their maximum values near the wall. Figure IS shows the profile of the resolv- 
able turbulent shear stress, <uv>, It can be seen that in the regions away from the walls the profile of 
<uv> does not follow the theoretical line, This indicates that the statistically stationary state has not 
been reached completely. Note that near the wall viscous stresses are Important, and the total shear stress 
must balance the gross pressure gradient. Moreover, in the present calculations, the subgrid-scale shear 
stresses are significant only very near the wall (y+ < 10). in Fig. 20, profiles of the intensities are com- 
pa red with seme of the available experimental data In the vicinity of the wall, .'he agreement of the computed 
<u" 2 >> < and <w 2 >*' 2 with the data Is satisfactory. However, as was also the case In Ref. 9, near the wall, 
a significant portion of <v 2 > 1 ' 2 seems to reside In subgrid-scale motions, This is consistent with our 
previous observation that vj is still excessively large near the walls. 

,, .Figure 21 shows_ the resolvable portions of the pressure velocity- gradient correlations, <p( ju/sx)n 
<jS(av/3y)>, and <p(aw/az)> In the vicinity of the wall (y + < 100, t » 2.0j. These terms are responsible for 
the exchange of energy between the throe components of resolvable turbulence kinetic energy, they are of 
particular Interest to turbulence modelers. Examination of these nroflles reveals that except In the Inroe- 
diate neighborhood of the wall ly' < 20), as expected, energy is I -ansferred from <u M2 » to <v 2 » end <w 2 >; 
that Is, <pUu/ax)> < 0, and *p(av/5y)>, -p(5w/»z)> > 0. On the other hand, as we approach the wall, a sig- 
nificantly different behavior can be noticed. Specifically, there is a relatively large rate of energy 
transfer from <v 2 >, whereas there Is a large energy transfer to <w 2 >. Tills rather unexpected result »s 
consistent nonetheless with our previous discussions of the fluid motions very close to the wall (Sec. 7), 

For example, Fig. 15a shows high-speed fluid approaching the wall and spreading laterally, resulting In 
relatively large energy transfer from <v 2 » to <i» 2 >. On the other hand, the momentum transfer from the 
lateral to the normal directions, which results in ejection Df fluid elements away from the wall, involves 
the nonenergetic (slow moving) fluid in the Immediate neighborhood of the wall. Thus, there Is a net energy 
transfer from <v-> to sw 2 *, as shown In Fig. 21. 

It should be mentioned that, In general, the values of the pressure velocity-gradient correlations 
computed In the present study are significantly higher than the earlier results using a much coarser grid 
(Ref. 9). This may indicate that a substantial portion of the pressure-strain correlation Is due to small- 
to-roadlum turbulence scales. To confirm this observation, several computations were carried out with differ- 
ent filter widths. The results of the calculations tend to support this observation. Thus, at present, 
and In the absence of a bettor subgrfd-scale turbulence theory, the computed pressure-strain correlations 
should be Interpreted qualitatively. It should be mentioned, however, that the large-scale flow structures 

presented in the previous section arc rather Insensitive (qualitatively) to the different filter widths and 

subgrid-scale models used. 

Before concluding this section, we turn our attention again to the subgrid-scale made! used 1r, the 
present study, to better resolve the relatively small turbulence scales in the vicinity of the walls, the 

present calculations were carried out for the case of a relatively low Reynolds number turbulent channel 

flow (Re T » 640, Re * 13,800), Therefore, the subgrid-scale turbulence Reynolds number defined in Sec. 3 is 
considered to be low in the regions away from the wall and very Tow In the vicinity of the walls. As was 
mentioned In Sec. 3, the arguments used In constructing this model arc valid only at a very high Reynolds 
number. Numerical results of McMillan and Farziger (Ref. 30) also show that Smagorinsky's model Is more 
appropriate at high Reynolds numbers. Thus, a low Reynolds number correction seems to be necessary. Note 
that because of the use of a much finer grfd In this simulation than that used In Ref. 9, the effective 
subgrid-scale turbulence Reynolds number Is lover than that In Ref. 9. In addition, because of the quasl- 
cycllc nature of turbulent channel flow (bursts, sweeps, etc.) the present calculations seem to indicate 
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that a subgrld-sealo model that has a bettor rojponso to tho tlroo history of the flow (a dynamic model) than 
the simple eddy viscosity model uiod here may bo necessary. TfiTs ft' necessary for a proper long-time Inte- 
gration of the governing equations. Integrating an additional equation for subprld-scale turbuTcnof energy 
Is an attractive possibility. In tho Interim, however, wo have round that selective filtering of the excess 
small-scale turbulence may oe adequate. 

9. CONCLUSIONS 

In this study, the three-dimensional time-dependent equations of motion have been numerically Integrated 
for the case of fully-developed turbulent channel flow. Tho calculations were carried out on tho ILLJAC IV 
computer with 64 mesh points In each of the spatial directions. Detailed flow patterns were studied by 
examining contour plots of typical Instantaneous velocity and vorticlty fields. In summary: 

1. The wall layer consisted of cohoront structures of Sow-speed and high-speed streaks alternating In 
the spanwise direction. These structures arc absent In the regions away from the wall. In addition, contour 
Plots of velocities In a typical y-z plane revealed tho existence of well -organized flow patterns In tho 
wall region, 

j. Hot spots, small localized regions of v^ry large values of turbulent shear stress, G'’v, were fre- 
quently observed. Very close to the wall, thosi hot sp.**s wore associated with 8" > 0 and 9 * 0 (sweep)*, 
away from tho wall, they were duo to fl" < 0 ana V > 0 (burst). In the central regions of tho channel, 
bursts from both halves of tho channel were the sources of the hot spots. 

3. No evidence of a direct relationship between streaks and stroamwisu vorticlty ’ x was observed 

In tho present simulation: very close to the wall, was not tho result of large-scale revolving fluid 
motions but was rattier due to the spanwise velocity gradient, (aw/ay), Though strong vortical regions were 
observed away from tho wall (y* - 30), attained Its maximum value at the wall, 

4. The profiles of the pressure velocity-gradient correlation showed a significant transfer of energy 
from the normal to tho spanwise component of turbulent kinetic energy In the Immediate neighborhood of tho 
wall (the "splattlng" effect), A large portion of the pressure-strain correlations appears to be due to 
small to medium scales of turbulent motions. 

me work presented here is still In progress and much more remains to be dono. In particular, a more 
refined model that depicts the dynamic nature of tho subgrid-scale motion may become necessary. Also, more 
mesh points, especially in the spanwise dfrectlon, are required in order to resolve the streaks at their 
proper scale. A computation with twice as many grid points as in tho present calculation (64 « 64 * 128) 
will be carried out In the near future. 

It Is hoped chat tills paper has demonstrated some of tho capabilities of LES as a research tool for 
Studying tho mechanics and structure of turbulent boundary layers. The authors believe that LES will make 
Important contributions to the study of turbulent flows by supplementing tho experimental data. 
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Fig. J. Contour plot of u"» In the «■: plane at y* • 16. 



Fig. 4. Contour plot of u“» In the x-z plane at y* ■ 90 
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1g. 6. Contour plot of u"& In the *-* plane at y/a • 0. 73 


Ftg. 6. Contours of the streanwlsc vortlclty tn the w plane at y* • 16. Note that the 
pattern' do not exhibit elongated structures In the x-dlrectlon. 
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OYI INJECTED AT THE MALL MILL BE COLLECTED MERE AND Fig. 17. Mean-velocity profile obtained with the 
LlFTEDurwARO sharp cutoff model (Ref 26). 


(•) Cross-sectional view of spanwlse velocity In 
y-l pUn#. 



fig. 18. Profile! ot horizontally averaged resolv- 
able turbulence Intensities. 
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Fig. 19. Vertical profile of horizontally averaged 
resolvable turbulent shear stress. 
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Flj. 70. Comparison of the horizontally averaged resolvable turbulence Intensities ullh experimental data. 



Fig. 71. Vertical profiles of horizontally averaqed resolvable pressure velocity qradlent correlations In 

the vicinity of the wall (y 4 * 100). 





